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ABSTRACT 


The  theoretical  basia  and  performance  characteristics  of  two  new  methods 
for  the  computation  of  the  coefficients  of  the  terms  of  asymptotic  expansions  at 
reentrant  comers  from  finite  element  solutions  are  presented.  The  methods,  called 
the  contour  integral  method  and  the  cutoff  function  method,  are  very  efficient: 
The  coefficients  converge  to  their  true  values  as  fast  as  the  strain  energy,  or  faster. 

In  order  to  make  the  presentation  as  simple  as  possible,  we  assume  that  the 
elastic  body  is  homogeneous  and  isotropic,  is  loaded  by  boimdary  tractions  only 
and,  in  the  neighborhood  of  the  reentrant  comer,  its  boundaries  are  stress  free. 
The  methods  described  herein  can  be  adapted  to  cases  without  such  restrictions. 
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1.  INTRODUCTION. 

Our  understanding  and  ability  to  control  the  error  of  approximation  in  finite 
element  computations  has  advanced  very  substantially  in  the  last  two  years.  Two 
developments  are  especially  important:  (l)  Under  assumptions  which  are  generally 
satisfied  in  engineering  computations,  we  are  now  able  to  realize  exponential  rates 
of  convergence  in  any  quantity  of  interest;  and  (2)  using  feedback  information  from 
finite  element  solutions  we  can  ensure  that  the  error  in  all  quantities  of  interest 
is  small.  In  this  paper  we  describe  two  new  methods  for  the  computation  of 
generalized  stress  intensity  factors  in  two  dimensional  problems  of  linear  elasticity 
and  demonstrate  their  performance  by  examples. 

We  will  restrict  our  attention  to  the  displacement  formulation.  We  denote 
the  solution  domain  by  0  and  the  displacement  vectors  defined  on  n  by  u  or 
(equivalently]  {u}: 

u  s  {u}  y)  u,(*,  y)}.  (1.1) 

We  denote  the  work  done  on  the  elastic  body  by  the  stresses  corresponding 
to  u  when  the  elastic  body  is  subjected  to  displacements  v  by  B(u,v);  the  strain 
energy  by  ll(u)  and  the  potential  energy  by  !!(»).  Note  that  2U(u)  =  B(u,u).  We 
denote  the  set  of  all  <Usplacement  vector  functions  defined  on  0  for  which  U(u)  <  oo 
by  E((l]  and  associate  the  energy  norm  with  this  set: 

ll«lls(o)  ='  (1.2) 

We  denote  the  set  of  all  admissible  displacement  functions  by  .^(Q).  By  definition, 
a  displacement  f\mction  u  is  admissible  if  it  has  finite  strain  energy  and  on  those 
boundary  segments  of  n  where  one  or  both  displacement  components  is  prescribed 
the  components  of  u  equal  the  prescribed  displacement  components  values.  That 
is: 

^(n) n* {u|ff€^(n):  u.{P)  =  u.(n  Pedni^^-,  u^(p)  =  u,(P),  p&an[°'>]  (1.3) 


where  u,,  Uy  are  prescribed  displacement  components  and  (resp. 

represents  those  boundary  segments  where  u,  (resp.  tiy)  is  prescribed.  The  exax:t 
solution  of  the  generalized  formulation  of  a  problem  of  plane  elasticity  Sex  satisfies: 

n{SBx)  -  nun  n{S).  (1.4) 

a€i{o) 

In  finite  element  analysis  we  seek  to  approximate  Ssx-  We  do  this  by  con¬ 
structing  a  finite  element  mesh  A  on  0.  Each  quadrilateral  (resp.  triangular)  finite 
element  is  mapped  onto  a  standard  quadrilateral  (resp.  triangular)  finite  element 
by  suitable  mapping  fimctions.  We  define  sets  of  basis  functions  on  the  standard 
triangular  and  quadrilateral  elements  so  that  any  polynomial  of  degree  p  or  less 
defined  as  the  standard  element  can  be  written  as  a  linear  combination  of  the  basis 
function. 

The  polynomial  basis  fimctions  defined  on  the  standard  element  are  mapped 
onto  the  elements  of  the  mesh  A  and  are  joined  to  form  a  set  of  basis  functions 
on  n.  These  basis  functions  are  continuous  across  interelement  boundaries  but  no 
restrictions  are  imposed  on  their  derivatives  across  interelement  boundaries.  Thus 
the  basis  functions  are  exactly  and  minimally  conforming.  The  basis  functions 
defined  on  0  are  characterized  by  the  mesh  A,  the  polynomial  degree  p  and  the 
mapping  fimctions  Q.  The  set  of  functions  that  can  be  expressed  by  linear  com¬ 
bination  of  the  basis  function  is  denoted  by  5^(0,  A,  Q)  and  called:  finite  element 
space.  Because  the  basis  functions  are  continuous  across  interelement  bound¬ 
aries,  5^(0,  A,  Q)  is  a  subset  of  £(0).  We  denote  the  set  of  admissible  functions  in 
5'’(n,  A,Q)  by  5'’(n,  A,Q).  The  number  of  basis  functions  in  S^{Cl,^,Q)  is  called  the 
number  of  degrees  of  freedo' 3  and  is  denoted  by  N.  The  finite  element  solution 
Sfb  satisfies: 

ff(arB)  —  nus  (1-5) 

«€S»(n,A,o) 

and  has  the  following  property: 


llwsx  -  =  min  l|»sx  -  ulls(n)-  (1-®) 

«€5mo.a,o) 

The  finite  element  method  selects  upg  from  5'’(n,A,Q)  on  the  basis  of  criterion 
(1.6).  We  can  reduce  the  error  of  approximation  by  mesh  refinement,  increase  of 
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the  polynomial  degree  of  elements,  or  a  combination  of  both.  These  are  called 
extension  processes.  If  the  polynomial  degree  of  elements  is  fixed  and  the  error  of 
approximation  is  reduced  by  mesh  refinement  then  the  process  is  called  h-extension 
(h  refers  to  the  size  of  the  elements).  If  the  mesh  is  fixed  and  the  polynomiad  degree 
of  elements  (p)  is  increased  then  the  process  is  called  p-extenaion.  We  remark  that 
neither  the  mesh  refinement  nor  the  distribution  of  polynomial  degrees  has  to  be 
uniform.  There  is  a  substantial  improvement  in  performance"'  and  no  significant 
increase  in  computational  overhead  if  properly  designed  meshes  are  used  instead 
of  uniform  meshes.  On  the  other  hand  there  is  no  significant  improvement  in 
performance  but  there  is  a  significant  increase  in  computational  overhead  if  graded 
p-distributions  are  used  instead  of  imiform  p>distributions.  For  this  reason  we  will 
use  \miform  p-distributions,  characterized  by  a  single  number,  p.  When  the  error 
of  approximation  is  reduced  so  that  mesh  refinement  is  combined  with  increasing 
p  then  the  extension  process  is  called  b~p  extension. 

In  the  case  of  h-  and  p-extensions  the  rate  of  convergence  is  algebraic.  That 
is; 


I  «SX  -  ttFS  (|s(0)  ^ 


(1.7a) 


where  N  is  the  number  of  degrees  of  freedom;  k  and  /S  are  positive  constants, 
independent  of  N  but  dependent  on  ubx  the  finite  element  meshes.  In  the 
case  of  h-p  extensions,  when  proper  mesh  refinement  is  used  in  conjunction  with 
p-extension,  the  rate  of  convergence  is  exponential: 


ISBX-Srs  IIk(o)  ^ 


(1.76) 


where  k,  7  and  ff  are  positive  constants.  Under  certain  assumptions  which  are 
generally  satisfied  in  engineering  applications  9  >  1/3. 

Although  it  is  difficult  to  write  computer  prograims  that  automatically  design 
meshes  and  assign  polynomial  degrees  so  that  the  rate  of  convergence  (1.7b)  is  re- 
zdized  for  all  AT,  very  similar  performance  can  be  achieved  when  properly  graded, 

*  In  this  context  ’‘performance’  meane:  rate  of  change  of  error  in  energy  norm  mth  respect 
to  the  number  of  degrees  of  freedom. 


fixed  meshes  are  iised  in  conjunction  with  p>-extension.  In  this  case  nearly  exponen¬ 
tial  convergence  rate  is  achieved  for  low  N  but  convergence  slows  to  an  algebraic 
rate  for  high  N.  Clearly,  the  mesh  should  be  designed  so  that  the  desired  level 
of  precision  is  achieved  before  the  rate  of  convergence  slows.  Properly  designed 
meshes  are  such  that  the  sizes  of  elements  decrease  in  geometric  progression  with 
a  common  factor  of  about  0.15  toward  points  of  stress  singularity.  Selection  of  the 
mesh  and  the  polynomial  degree  of  elements  depends  on  the  accuracy  one  wishes 
to  achieve.  For  details  and  examples  see  [1-8]. 

We  are  now  in  position  to  outline  the  central  point  of  this  paper:  In  [9]  a 
new  approach  for  the  computation  of  ftmctionals  from  finite  element  solutions 
was  presented.  In  this  approach  a  desired  functional  value  fiUFs)  (e.g.  a  stress 
component  at  a  point;  stress  intensity  factor,  etc.)  is  computed  from  an  expression 
which  is  of  the  form: 

#(«rs)  =  B{ufs,  v)  (1.8) 

where  v  is  a  suitably  chosen  fxmction,  called  extraction  function.  Methods  of  this 
type  are  called  extraction  methods.  In  [10]  it  was  shown  that  the  error  in  the 
fimctional  value  computed  by  extraction  methods  can  be  written  as: 

l^(«Bx)  -  ^(«ps)|  <  II^BX  -  WBBllBCnjjjwBX  -  WBb|Ib(o)  (1-9) 

where  (3fe  is  the  finite  element  solution  of  an  auxiliary  problem,  the  exact  solution 
of  which  is  ubx-  In  the  auxiliary  problem  the  domain  and  the  mesh  are  the  same  as 
in  the  original  problem  but  the  loading  is  computed  from  the  extraction  fimction 
V.  The  auxiliary  solution  serves  theoretical  purposes  only:  It  is  not  necessary 
to  know  uex  or  ufb  in  order  to  use  (1.8)  for  computing  ^iuBs).  It  is  possible 
to  select  the  extraction  ftmction  «  so  that  \\uex  -  "bbIIbco)  0  not  slower  than 
||«BX  -«FB||B(n)  -*0as  N  ooinan  orderly  extension  process.  Therefore  functional 
values  computed  by  extraM:tion  methods  have  the  szune  order  of  accuracy  as  the 
strain  energy  or  better.  Functional  values  computed  by  extraction  methods  exhibit 
superconvergence  [10,11]. 

In  this  paper  we  describe  procedures  for  the  extraction  of  the  amplitude  of 
stress  singular  terms  associated  with  reentrant  comers  in  plane  elasticity.  We  as¬ 
sume  that  the  plane  elastic  body  is  homogeneous  and  isotropic;  loaded  by  bound¬ 
ary  tractions  only;  in  the  neighborhood  of  the  reentrant  comer  the  boundaries  eire 


stress  free,  zuid  in  (1.3)  u,  =  0,  u„  =  0,  however  the  method  described  herein  can 
be  generalized  to  other  problems  such  as  fixed-free  boundary  conditions;  inter¬ 
section  of  material  interfaces  with  external  boundaries;  etc.  We  demonstrate  the 
performance  of  the  extraction  procedures  on  the  basis  of  three  test  problems. 

2.  THE  EXACT  SOLUTION  IN  THE  NEIGHBORHOOD  OF 
REENTRANT  CORNERS. 

In  the  neighborhood  of  reentrant  comers  (Fig.  2.1)  we  examine  Airy  stress 
functions  of  the  form: 

U  =  U{r,e)  =  r>^*^F{e).  (2.1) 


Fig.  2.1.  Reentrant  comer.  Notation. 


Because  U  satisfies  the  biharmonic  equation: 

f  1  ^  (£_  ld_  £_\ 

^  r  dr  ^  dd^ )  ^  r  dr  ^  d9^  ) 


(2.2) 


F{6)  must  satisfy: 

F""  +  2  (A®  +  1)  F"  +  (A»  -  1)»  F  =  0  (2.3) 

where  the  primes  represent  differentiation  with  respect  to  9.  The  general  solution 
of  (2.3)  for  A  /  0  and  A  #  ±  1  is: 


F(9)  =  ai  cos(A  —  1)9  +  aa  co8(A  +  l)^  +  03  8iii(A  —  l)tf  +  04  8in(A  +  l)fl. 


(2.4) 


We  need  to  detennine  A  and  combinations  of  (t  =  i,  2,  3,  4),  such  that  the  edges 

that  meet  at  the  reentrant  corners  are  stress  free.  From  the  stress  function  U  we 
have: 


(2.5a) 


and: 


la^U  idU  , 


(2.56) 


Along  the  reentrant  edges  (at  6  =  ±a/2)  we  have  a*  =  =  o.  Consequently,  from 

(2.4)  and  (2.5a, b),  after  straightforward  algebraic  manipulation,  we  have: 


C08(A-1)| 

co8(A  +  1)| 

0 

0 

’  ax  ' 

A»m(A-  1)| 

8in(A  +  1)| 

0 

0 

<*2 

0 

0 

8iii(A-  1)| 

8in(A  +  1)| 

< 

Os 

0 

0 

C08(A-1)| 

A  co8{A  +  t)^ 

■04  , 

>  =0 


(2.6) 


where: 


A  = 


1- A 


1  A 


(2.7) 


Note  that  ai,  are  independent  of  as,  a4.  A  nontrivial  solution  exists  only  if  the 
corresponding  determinants  vanish: 


coe(A  -  l)^8in(A  +  1)^  -  A8in(A  -  l)^co8(A  +  1)^  =  0 
«  *  *  2 


8in(A  -  1)  ^  cob(A  +  1)  ^  —  a  co8(A  —  1)  ^  8m(A  +  1)  —  =  0 
2  2  2  2 


(2.8a) 

(2.86) 


which  can  be  simplified  to: 


8inAa  +  A8Uia=:0 
8iii  A  a  —  A  8ia  a  =  0. 


(2.9a) 

(2.96) 


See  also  [12-16].  Observe  that  if  A  is  a  solution  of  (2.9a)  or  (2.9b),  then  -A  is  also 
a  solution,  however  the  corresponding  stress  field  has  finite  strain  energy  only  if 
the  real  pairt  of  A  is  greater  than  zero. 


J. 


Assume  now  that  Aj*'  (»  =  1,2,...)  is  a  solution  of  (2.9a)  and  A^’  is  real  and 
simple.  Then,  from  (2.6): 


ai  coe(Al‘'  -  1)-  +  aj  co8(Aj^’  +  1)  -  ^  0 


ai  8in(Af  ’  -  l)y  +  oa  co«(Ap’  +  1)|  =  0 


(2.10a) 


(2.106) 


where:  = 


l-A* 


1+Ar 


Let  us  define: 


^  ^  C08(a|‘‘-1)-  A,'^' 8m(A)^’ - 1)- 

•  «i“  c08(A<*>  +  1)|  sm(A<^>  +  l)| 


With  this  notation  (2.1)  can  be  written  as: 


U  =  (co8(A.*‘‘  -  1)  co8(a'‘>  +  1)6). 


(2.12a) 


Similarly,  if  a|^'  is  a  real,  simple  root  of  (2.9b),  then  the  corresponding  stress 


function  is: 


where: 


U  =  (8in(Ap>  -  1)  ^  ^  sm(A|’^  +  l)  o) 


(2.126) 


8m(Af  *  -  1)|  ^  co8(a|^^  -  1)| 

8in(A|’>  +  1)|  1  +  a|’>  coe(Ap’  +  1)| 


Note  that  (2.12a)  is  symmetric  with  respect  to  9,  whereas  (2.12b)  is  antisym¬ 
metric.  From  (2.12a,b)  expressions  from  stress  and  displacement  fields  can  be  de¬ 
rived.  These  expressions  are  most  conveniently  obtained  by  first  writing  (2.12a,b) 
in  terms  of  the  complex  variable  2,  and  then  iising  the  method  of  Muskhelishvili 
[15].  The  stress  and  displacement  fields  corresponding  to  (2.12a)  are  called  Mode 
1  fields.  The  stress  and  displacement  fields  corresponding  to  (2.12b)  are  called 
Mode  2  fields.  Specifically,  the  Mode  1  displacement  components,  up  to  rigid 
body  displacement  and  rotation  terms,  are: 

[(*-«l‘’(A|'’  +  l))co8A'‘>(?-A|'>co8(Ai^'  -2)«]  (2.14a) 

=  [(*  +  <^i‘’(Ar’  +  l))8mAp>(»  +  A|‘>8m(Al'>-2)tf]  (2.146) 
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which  can  be  written  in  the  form: 


(2.14c) 


In  (2.14a,b,c)  G  is  the  modulus  of  rigidity,  and  k  depends  only  on  Poisson’s  ratio. 
For  plane  strain: 

<c  =  3-4i/  (2.15o) 

and  for  plane  stress: 

Z-v 


K  = 


1  +  t/ 


(2.156) 


The  Mode  1  stress  tensor  components  are: 


i(2-<3f>(Aj‘)  +  l)] 

1 

11 

(2  +  <?r’(Ai'>  +  l)) 

'*v*  ''t  '  1 

(A<‘>-l)8iii(Al^>- 

The  Mode  2  displacement  components,  up  to  rigid  body  displacement  and  rotation 
terms,  are: 

^  [(«  -  +  1))  sin(A.I^>  -  2)  $]  (2. 17o) 

[('*  +  1))  -  2) »]  (2.176) 

which  can  be  written  in  the  form: 


(2.17c) 


The  Mode  2  stress  tensor  components  are: 

<^iV  =  a]’’  [(2  -  <?|^’(Af  ’  +  1))  sui(Af  >  -  1)  ^  -  (Ap’  -  l)am(Ap’  -  3)  <»](2.18a) 

=  aI^’  [(2  +  <3|’’(A*.’’  +  1))  sin(A|^’  -  1)  ^  +  (aJ“>  -  l)8m(A|®'  -  3)  «]  (2.186) 

=  -Aj"’’  [(Aj^’  -  1)  co8(Ai^>  -  3)  S  +  Q|’’(a|’>  +  1)  co8(Aj®'  -!)«].  (2.18c) 

In  fracture  mechanics  we  have:  a  =  2ir.  In  this  case  (2.9a,b)  are  identical 
(sm2XT  =  o),  and  all  roots  are  real  and  simple: 

Al‘'=Ai"  =  4  4  ±2,  4... 
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(2.19) 


Assuming  traction  free  crack  surfaces,  finite  strain  energy,  and  neglecting  rigid 
body  displacement  2uid  rotation  terms,  in  the  neighborhood  of  the  crack  tip  any 
solution  can  be  written  in  the  form: 


(2.20) 


This  infinite  series  converges  absolutely  for  r  <  ro  for  some  ro  >  0.  The  coefficients 
and  are  called  genenlixed  stress  intensity  factors.  The  coefficients 
and  A^^^  2n^e  related  to  the  Mode  1  and  Mode  2  stress  intensity  factors,  Ki  and 
Kii,  as  follows: 


(2.21) 


Methods  for  the  determination  of  the  generalized  stress  intensity  factors  are  dis¬ 
cussed  in  the  following  sections. 

When  a  /  2ir,  then  not  ail  roots  are  real,  and  mtiltiple  roots  (both  real  and 
complex)  can  exist.  The  cases  of  complex  roots  and  multiple  real  roots  are  beyond 
the  scope  of  this  discussion.  From  the  point  of  view  of  engineering  analysis,  the 
solid  angles  of  270,  240,  225  and  210  degrees  are  of  the  greatest  importance  (in 
addition  to  the  360  degrees  already  discussed).  In  these  cases  the  smallest  roots 
are  real.  Their  values  are  shown  in  Table  2.1. 


Table  2.1.  Smallest  positive  roots  of  (9a,b) 
for  selected  solid  angles  (a). 


a 

A‘^) 

A‘») 

360* 

0.500000 

0.500000 

270* 

0.544484 

0.908529 

240* 

0.615731 

1.148913 

225* 

0.673583 

1.302086 

210* 

0.751975 

1.485812 

When  the  solid  angle  is  greater  them  257.45  degrees,  then  >  i  and  the 
Mode  2  stress  components  are  botmded.  The  Mode  1  stress  components  are  un¬ 
bounded  when  a  >  180*. 


I 

3.  BETTI’S  LAW.  THE  PATH  INDEPENDENT  INTEGRAL  /r.(«,t7). 


We  denote  the  strain  tensor  components,  corresponding  to  the  displacement 
field  u(z,y)  by  £1"^  The  strain-displacement  relations  are: 


c 


dus  9m, 
dy  9i 


(3.1) 


Similarly,  we  denote  the  stress  tensor  components  corresponding  to  displacement 
field  u(i,  y)  by  and  rlj* .  The  stress-strain  relations  are: 


=  {£<-)} 


(3.2) 


where  {«<“>}  =  and  [E]  is  a  symmetric, 

positive  definite  matrix,  called  the  material  stiffness  matrix.  We  have  assumed 
that  the  initial  strain  is  zero.  We  shall  assume  also  that  the  body  forces  are  zero 
and  that  u  is  such  that  the  coresponding  stresses,  computed  from  (3.1)  and  (3.2), 
satisfy  the  equilibrium  equations  with  zero  body  forces: 


£^  +  £j!L=0; 

dz  dy 


dz  dy 


(3.3) 


We  denote  the  direction  cosines  of  the  normal  to  the  boundary  of  the  plane 
elastic  body  at  boimdary  point  P  by  n„  n,.  The  traction  vector  components  in 
terms  of  the  stress  components  at  point  P  are: 

+  (3.4a) 

i;(“)  =rWn.  +  (3.4i) 


Let  «  =  {«}  =  {wa(z,y)  w,(z,y)}  be  an  wbitrary  displacement  vector  field  and 
assume  that  the  strain  energy  associated  with  v  is  finite  on  n.  Assxime  further  that 
tractions  are  specified  along  the  entire  boundary  of  0.  We  denote  the  boimdary 
of  n  by  dO.  In  the  absence  of  body  forces,  thermal  loading  and  elastic  constraints 
the  principle  of  virtual  work,  states  that: 

jj^{ti'’Y{ai-)}dxdy  =  y'^{€<*)}*’[£l{€(-)}«izdy  =  (3-5) 


I 

I 
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holds  for  any  v.  v  is  called  virtual  displacement.  The  strain  components  are 
related  to  v  by  (3.1),  Because  {E\  is  symmetric,  (3.5)  can  be  written  as: 

11^  +2;‘v„)da  (3.6) 

which  is  the  same  as: 


+  r<*) 

-r 


(3.7) 


Applying  Green’s  lemma  and  using  (3.4a, b)  we  have: 


/  (Ti«>ti.  +  Ii*>u,)da- /  (ri“)t».  +  i;(“)t,,)ti5  = 
Jan  Jan 


When  the  stresses  corresponding  to  both  u  and  v  satisfy  the  equilibrium  equations 
with  the  body  forces  equal  to  zero,  that  is  (3.3)  then  we  have  Betties  law: 

f  +  Tj*)tv)  ds=  f  (!;<“)«,  +  T(»)t,^) da  (3.9a) 

Jan  Jan 

where  the  integration  is  countreclockwise  around  n.  We  denote  the  normal  and 
tangential  traction  vector  components  respectively  by  r„  and  Tt  and  the  normal  and 
tangential  displacement  vector  components  by  un  and  ut.  Because: 

etc,,  Betti’s  law  can  be  also  written  in  the  following  form: 

/  (Ti*) Un  +  T, <•’«,)  (ia  =  /  (7l“)v„  +  r/’‘>u.)<f3.  (3.96) 

Jan  Jan 


Let  us  now  consider  a  subdomain  of  n,  denoted  by  0*,  in  the  neighborhood 
of  the  reentrant  comer,  n*  is  bounded  by  two  continuous  curves,  and  F,,  and 
the  reentrant  edges,  as  shown  in  Fig.  3.1.  Let  {u}  and  {v}  be  two  displacement 
fields,  both  satisfying  the  equilibrium  equations  with  body  forces  and  initial  strains 
equal  to  zero,  and  the  stress  free  boundary  conditions  along  the  reentrant  edges 
(ffe  =  r,*  =  0  at  =  ±a/2).  Then,  according  to  Betti’s  law: 
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i*m  t.  * 


.**U 


Fig.  3.1.  The  subdomain  n*. 

f  {Ti:^Un  +  Ti^K)da+  f  (Ti^^un +  Ti'^nt)ds  = 

Jr\  Jr; 

f  +  Ti^^iH)ds.  (3.10) 

Jrx  Jri 

We  integrate  arotmd  0*  in  the  counterclockwise  direction.  Therefore  integration 
along  r;  is  in  the  cotinterclockwise  direction,  while  integration  along  rj  b  in  the 
clockwbe  direction  with  respect  to  the  reentrant  comer.  Let  us  reverse  the  sense 
of  the  direction  of  the  integral  along  H  so  that  both  contour  integrab  are  in  the 
cotmterclockwbe  direction  with  respect  to  the  reentrant  comer.  We  will  use  the 
symbol  f;  to  indicate  integration  along  r;  in  the  counterclockwbe  direction  with 
respect  to  the  reentrant  comer.  In  thb  case  we  have: 

Jfi  Jfi 

f  +  /  (4*)u„  +  r/’>u.)«i3.  (3.11) 

^r;  Jri 

Therefore,  the  integral: 

/r.(t?,5)=/  (3l->w„  +  I^“‘w)*- /  (2;i*‘u„  +  3;f’’u.)* 

Jr*  Jr* 


(3.12) 


is  independent  of  the  path  P*.  Of  course,  P*  must  begin  on  the  reentrant  edge  and 
terminate  on  the  other  reentrant  edge,  as  shown  in  Fig.  3.1,  and  the  integration 
must  be  counterclockwise  with  respect  to  the  reentrant  comer. 

When  the  stresses  corresponding  to  {»}  do  not  satisfy  the  equilibrium  equa¬ 
tions  and  the  stress  free  boundary  conditions  along  the  reentrant  edges,  then 
instead  of  (3.11)  we  have: 


/  {TirKn  +  Ti^^vt)d,-  f  +  = 

Jt\  Jtl 

f  (Tir^v^-^Ti^U)ds-  f  +  f  {T{:^u^  +  Ti•'>iH)ds- 

JT\  Jti 


The  sense  of  integration  along  P3  and  P4  is  shown  in  Fig.  3.1. 


4.  EXTRACTION  OF  STRESS  INTENSITY  FACTORS. 

In  the  following  we  consider  only  real  and  simple  roots  of  (2.9a)  and  (2.9b). 
This  is  the  case  in  linear  elastic  fracture  mechanics.  We  present  two  algorithmic 
procedures  for  the  computation  of  coefficients  (t  =  1,  2,...;  m  =:  l,  2)  in  (2.20). 
The  restriction  to  real  and  simple  root  is  not  essential  however,  and  all  coefficients 
of  asymptotic  expansions,  similar  to  (2.20),  can  be  determined  by  the  methods 
described  in  this  section. 


4.1.  The  contour  integral  method. 


Let  be  a  circle  of  radius  p  centered  on  the  crack  tip,  and  assume  that  p  is 
sufficiently  close  to  the  crack  tip  so  that  the  exact  solution  {usx}  is  represented 
by  (2.20)  on  The  traction  vector,  corresponding  to  {usx},  can  be  written  as: 


•si  msl 


(4.1) 


where  >  0  is  a  real  and  simple  root  of  (2.9a)  and  a|‘^  >  0  Is  a  real  and  simple 
root  of  (2.9b). 


Let  -Ay  ’  be  a  negative  root  of  (2.9a)  or  (2.9b),  and  denote  the  corresponding 
displacement  function  by  {wL"*  }•' 

{«!.'*;}  =  (4.2) 

where  is  a  constant,  to  be  determined  later.  Note  that  (vi”^}  does  not  have 
finite  strain  energy  on  0  and,  therefore,  is  not  an  admissible  displacement  fimction 
on  n.  It  is  admissible,  however,  on  n*  (Fig.  3.1).  The  traction  vector  on  r^, 
corresponding  to  is: 

(ri/'^’’}  =  a)'**  (tf)}.  (4.3) 

We  now  show  that: 


f.4L7‘A<->c<"‘>(a)  ifi=ji 
1 0  otherwL 


=  j  and  m  =  I 


1 0  otherwise 

where  Tp  is  a  circular  arc  of  radius  p  centered  on  the  reentrant  comer,  cj^^la) 

depends  only  on  aI"*’.  The  proof  is  straightforward.  By  direct  evaluation  from  the 

« 

absolutely  convergent  series  (2.20),  we  have: 

=  (-=1.2)  («.) 


where: 


1  /  \ 

+  {»!'"* (»)r{T<:;(^)})  de.  (4.56) 

Because  /r,(«sxi  vL"’)  is  independent  of  Tp,  it  is  independent  of  p.  If  aJ”*  /  Aj"’,  then 
this  is  possible  only  if  the  integral  expression  (4.5b)  is  zero.  The  other  possibility 
is  that  aI*"^  =  Aj.”*,  but  n  jt  m.  In  this  case,  however,  the  integrand  in  (4.5b)  is 
the  sum  of  the  dot  products  of  symmetric  auid  antisymmetric  vector  functions 


therefore  the  integral  is  zero.  Denoting:  cl" '(a)  =  C','”;'"’,  we  have  (4.4)  and  letting 
we  have  from  (3.10): 

=  i,  ds  - ds  (4.6) 

where  Tj  is,  of  course,  arbitrary.  The  function  {«L7')  “  an  extraction  function 
for  Of  course,  we  do  not  know  {aax}  smd  {7^““)}.  We  therefore  substitute 
{ufc}  and  in  (4.6)  for  {ugx}  and  {T^****}  to  obtain  an  approximation  to 

vll”*).  If  r*  includes  an  external  boundary  where  the  imposed  tractions  are  known, 
then  the  imposed  tractions  can  be  used  in  (4.6).  Implementation  would  pose 
some  difficulties,  however,  and  computational  experience  has  shown  that  the  use 
of  instead  of  the  imposed  tractions  generally  yields  satisfactory  results. 

Examples  are  presented  in  Section  5. 

4.2.  The  cutoff  ftmction  method. 

Let  us  now  assume  that  both  contours  FI  and  P;  are  circular  arcs  with  radii 
Pi  and  Pa,  respectively,  Pi  <  pa.  We  define  the  following  extraction  function  for 

{wLf}«^(r){vL7’}  (4.7) 

where  {vi7^}  is  defined  by  (4.2),  and  ^(r)  is  called  the  cutoff  function  and  is  defined 


1  r  <  Pi 

^(r)=|  +  Pi<r<pa  (4-8) 

\P2-Pl/  \P2-Pl/ 

.0  r>p2. 

In  this  case  the  extraction  function  satisfies  neither  the  stress  free  boundary  con¬ 
ditions  on  the  reentrant  edges  nor  the  equilibrium  equations.  Therefore,  (3.13) 
must  be  used  for  the  computation  of  instead  of  (3.11).  Specifically,  if  we 
select  aL7^  as  before,  i.e:  a17*  =  where  cj”*'(a)  is  computed  from  (4.5b) 

with  {wL7*}  replaced  by  {t»L7*}.  we  have: 


T*[- 


dxdy. 


The  contour  integrals  must  be  evaluated  along  the  reentrant  edges  in  the  counter¬ 
clockwise  direction  with  respect  to  0*  (see  Fig.  3.1).  The  advantage  of  the  cutoff 
function  method  is  that  stresses  and  tractions  corresponding  to  {urg}  do  not  have 
to  be  considered.  For  this  reason  the  cutoff  function  method  is  more  accurate  than 
the  contour  integral  method. 

5.  EXAMPLES. 

The  contour  integral  and  cutoff  function  methods  can  be  implemented  in  any 
finite  element  computer  program.  However,  the  error  in  the  computed  data  is 
closely  related  to  the  error  in  energy  norm  and  therefore  the  design  of  the  finite 
element  space  A,Q)  [10,11].  Also,  the  cutoff  fimction  method  requires  the  use 
of  smooth  mapping  techniques.  The  following  example  problems  were  solved  by 
means  of  a  new  computer  program,  called  PROBE.  PROBE  implements  the  p-versjon 
of  the  finite  element  method,  that  is  sequences  of  finite  element  spaces  5>’(n,A,Q) 
can  be  conveniently  created  by  letting  p  =  l,2,...8  while  keeping  the  mesh  A  and 
the  mappings  Q  fixed.  Of  course,  the  mesh  and  the  mappings  can  be  changed 
also.  The  mapping  of  finite  elements  is  by  the  linear  blending  function  method. 
Curved  boundaries,  such  as  circles,  and  other  conical  sections  are  represented 
exactly  in  the  computation  of  stiffness  matrices  and  load  vectors.  PROBE  also 
permits  specification  of  loading  by  FORTRAiV-like  expressions,  hence  boundary 
tractions  specified  herein  were  represented  exactly  in  the  load  vector  computations. 
Twelve  Gauss  points  used  in  the  computation  of  load  vectors  independently  of  the 
polynomial  degree.  In  the  case  of  extraction  methods  12  Gauss  points  are  used  for 
evaluating  the  contour  integrals  in  and  144  Gauss  points  are  used  for  evaluating 
the  area  integrals,  independently  of  the  polynomial  degree.  Additional  details 
concerning  PROBE  are  available  in  [17|. 

Other  examples  of  the  application  of  extrau:tion  methods,  based  on  experi¬ 
mental  (research)  implementations,  are  available  in  [18,19]. 

5.1.  L-shaped  plane  elastic  body. 

We  consider  the  L-shaped  plane  elastic  body  of  thickness  t,  shown  in  Fig. 

5.1,  loaded  by  tr£M:tion8  corresponding  to  the  first  symmetric  and  antisymmetric 


a 


eigenfunctions  of  the  asymptotic  expansion  of  ubx  about  the  reentrant  comer. 
This  example  is  representative  of  plate  eind  shell  intersections  and  reentrant  corner 
problems  in  general. 


Fig.  5.1.  L-shaped  plane  elastic  body.  Mesh  design  (18  elements). 


We  select  u  =  0.3  and  assume  plane  striuji  conditions.  Therefore  k  s:  i.s; 
X[V  =  0.54448374;  ~  0.543075579;  =  0.90852919;  *  -0.218923236  [see  (2.15a), 

(2.9a, b),  (2.11)  and  (2.13)].  We  superimpose  the  tractions  corresponding  to  the 
stress  components  (2.16)  and  (2.18)  along  the  boimdaries  of  the  L-shaped  plane 
elastic  body  with  the  stress  intensity  factors  selected  so  that  *  A,  (A  is  arbi¬ 
trary)  and: 


(5.1) 


where  a  is  the  dimension  shown  in  Fig.  5.1.  In  this  way  the  exact  strain  energy 
can  be  computed: 


Z/(ubx)  =  6.77776914 


E 


(5.2) 


where  E  is  the  modulus  of  elasticity. 

In  the  contour  integral  method  (CIM)  the  integration  was  performed  along 
the  circle  of  radius  0.15  a  and  the  finite  element  solutions  were  computed  from  the 
elements  inside  the  circle. 

The  number  of  degrees  of  freedom,  the  computed  values  of  the  normalized 
strain  energy  and  the  normalized  stress  intensity  factors,  defined  by: 


are  listed  in  Table  5.1.  Of  course,  has  to  converge  to  1  and,  in  view  of  (5.3), 
a[*’  has  to  converge  to  2. 


Table  5.1.  L-shaped  domain.  Strain  energy  and 
the  normalized  stress  intensity  factors  a[^\  a[^*  computed  by 
the  contour  integral  method  (CIM)  aind  the  cutoff  function  method  (CFM). 


N 

li{uFB)E 

■*1 

1(3) 

1(3) 

■*1 

p 

(Aa*i  *)^t 

(CIM) 

(CFM) 

(CIM) 

(CFM) 

1 

41 

8A2072796 

1.18041 

0.95268 

2.43474 

2.29075 

2 

119 

6.74137580 

0.95418 

1.02177 

2.01352 

2.08422 

3 

209 

6.77029847 

1.02786 

1.00250 

2.02597 

2.02239 

4 

335 

6.77575144 

0.99014 

1.00073 

1.99801 

2.00437 

5 

497 

6.77683967 

1.00444 

0.99991 

2.00265 

2.00097 

6 

695 

6.77719530 

0.99784 

0.99985 

1.99939 

2.00022 

7 

929 

6.77736281 

1.00074 

0.99987 

2.00036 

2.00005 

8 

1199 

6.77746228 

0.99952 

0.99990 

1.99988 

2.00001 

oo 

oo 

6.77776914 

1.00000 

1.00000 

2.00000 

2.0000 

We  see  from  Table  5.1  that  the  stress  intensity  factors  computed  by  both 
methods  converge  strongly  and  obviously,  although  not  monotonically.  Greater 
accuracy  and  more  nearly  monotonic  convergence  is  exhibited  by  the  cutoff  func¬ 
tion  method  than  the  contour  integral  method,  nevertheless  both  methods  yield 
solutions  which  are  within  the  range  of  precision  normally  needed  in  engineering 
computations  at  p  =  2  or  p  =  3. 

We  have  plotted  the  relative  error  in  strain  energy  and  the  absolute  value  of 
the  relative  error  in  the  Mode  1  and  Mode  2  stress  intensity  faictors  computed  by 
the  cutoff  function  method  against  the  number  of  degrees  of  freedom  on  log-log 
scale  in  Fig.  5.2a  and  on  a  semilog  scale  (the  logarithms  of  the  relative  errors 
vs.  in  Fig.  5.2b.  These  choices  of  scale  are  motivated  by  estimates  (1.7a) 

and  (1.7b),  respectively.  These  diagrams  indicate  that  the  error  in  strain  energy 
behaves  differently  in  the  range  of  low  p  values  than  in  the  range  of  high  p  values. 
In  the  range  of  low  p  values  it  curves  downward  in  Fig.  5.2a  and  very  nearly 
follows  a  straight  line  path  in  Fig.  5.2b,  suggesting  that  estimate  (l.7b)  holds.  We 
will  refer  to  this  as  Phase  1.  In  the  range  of  high  p  values  it  follows  a  straight  line 
path  in  Fig.  5.2a  but  has  a  positive,  decreasing  curvature  in  Fig.  5.2b.  We  will 
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refer  to  this  as  Phase  2.  Transition  from  Phase  1  to  Phase  2  occurs  at  about  p  =  4. 
Phase  2  is  the  asymptotic  convergence  for  p-extensions:  Estimate  (1.7a)  holds 
with  p  =  in  this  case  p  =  1.089.  Phase  1  is  characteristic  for  h-p  extensions. 

Obviously,  the  rate  of  decrease  of  the  error  is  much  faster  in  Phase  1  than  in 
Phase  2.  We  can  extend  Phase  1  by  refining  the  mesh  further  so  that  the  elements 
are  graded  in  geometric  progression  toward  the  reentrant  comer  with  a  common 
factor  of  about  0.15.  Thus  the  next  layer  of  elements  around  the  reentrant  corner 
would  have  the  size  of  O.lS^a.  In  general,  the  mesh  should  be  designed  so  that 
the  desired  level  of  precision  is  achieved  in  Phase  1.  In  this  exsunple  the  relative 

errors  are  less  than  one  percent  in  Phase  1,  therefore  no  further  refinement  was 
necessary. 

The  convergence  path  of  the  Mode  1  stress  intensity  factor  follows  closely 
that  of  the  strain  energy,  whereas  the  the  Mode  2  stress  intensity  factor  converges 
faster.  This  point  is  discussed  in  [10,11]. 

5.2.  Edge  cracked  panel  problem  1. 

Let  us  now  consider  the  edge  cracked  panel  shown  in  Fig.  5.3.  We  assume 
plane  strain  conditions  and  Poisson’s  ratio  of  0.3.  In  this  case  =;  a(^*  =  i/2, 
therefore  Q*/*  =  -1/3  and  =  i.  Once  again  we  denote  the  thickness  of  the  panel 
by  t.  We  apply  tractions  along  sides  A,  B,  C,  D,  E,  F  of  the  edge  cracked  panel 
shown  in  Fig.  5.3  that  exactly  correspond  to  the  stresses  of  Mode  1  and  Mode  2 
stress  fields,  in  this  case,  using  appropriate  trigonometric  identities,  (2.16a,b,c) 

and  (2.18a,b,c)  can  be  written  in  the  following  form  (see,  for  example,  [20]): 
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where  -r  <  <  r.  The  Mode  2  stress  components  are: 
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The  generalized  stress  intensity  factors  and  are  related  to  the  stress 
intensity  factors  by  (2.21).  We  select  Kj  =  Ku  =  y/^ A  [A  is  arbitrary)  and  define 
the  normalized  stress  intensity  factors  and  as  follows: 


(5.6) 


In  this  way  the  computed  values  of  both  and  A^^^  converge  to  1  and  therefore 
it  is  easy  to  monitor  convergence  of  the  stress  intensity  factors.  The  exact  strain 
energy  is  known: 


Fig.  5.3.  Edge  cracked  panel  problem  1.  Mesh  design  (24  elements). 

The  number  of  degrees  of  freedom,  the  computed  values  of  the  normalized 
strain  energy  and  the  normalized  stress  intensity  factors  defined  in  (5.6)  are  shown 
in  Table  5.2. 

We  see  from  Table  5.2  that  again  the  cutoff  function  method  is  somewhat 
more  accurate  than  the  the  contour  integral  method.  In  the  case  of  the  cutoff 
fxmction  method  the  relative  error  falls  below  1  percent  at  p=3;  in  the  case  of  the 
contour  integral  method  the  relative  error  drops  below  1  percent  at  p=5. 

The  Relative  errors  in  strain  energy  juid  the  absolute  value  of  the  relative 
error  in  the  Mode  1  and  Mode  2  stress  intensity  factors,  computed  by  the  cutoff 


Table  5.2.  Values  of  and  computed  by  the 
contour  integral  and  cutoff  function  methods. 


p 

N 

U(iiFB)E 

at 

(CIM) 

(CFM) 

(CIM) 

(CFM) 

1 

53 

9.5959919 

1.16972 

0.91371 

1.31197 

1.00163 

2 

155 

10.4114292 

0.92515 

1.02702 

0.92840 

1.04058 

3 

273 

10.5085887 

1.03787 

1.00277 

1.04413 

1.00637 

4 

439 

10.5302827 

0.98483 

1.00025 

0.98360 

1.00134 

5 

653 

10.5354556 

1.00609 

0.99943 

1.00737 

0.99992 

6 

915 

10.5374132 

0.99650 

0.99944 

0.99639 

0.99978 

7 

1225 

10.5384252 

1.00083 

0.99954 

1.00125 

0.99978 

8 

1583 

10.5390566 

0.99907 

0.99964 

0.99917 

0.99982 

oo 

oo 

10.5412281 

1.00000 

1.00000 

1.00000 

1.00000 

function  method,  are  plotted  against  the  number  of  degrees  of  freedom  on  log-log 
scale  in  Fig.  5.4.  The  inverted  S-curve,  typical  of  p-convergence  when  strongly 
graded  meshes  are  used,  is  clearly  vbible.  The  stress  intensity  factors  converge  at 
about  the  same  rate  as  the  strain  energy.  The  curves  appear  to  enter  Phase  2  at 
p  =  7  or  p  =  8.  The  results  are  in  agreement  with  the  theoretical  estimate  given  in 
[10]. 

POLYNOMIAL  DEGREE 
I  2  3  4  5  6  7  8 


Fig.  5.4.  Convergence  of  the  strain  energy  and  the  Mode  1  and  Mode  2 
stress  intensity  factors  computed  by  the  cutoff  function  method. 


Finally,  let  us  examine  the  vector  length  of  stress  intensity  factors  which 


can  be  computed  by  the  energy  release  rate  method  also  known  as  the  stiffness 
derivative  method.  Because  we  have  selected  for  this  test  problem  Ki  =  Ku  = 
V^i4,  we  define  the  computed  values  of  of  the  normalized  vector  length  of  the 
stress  intensity  factors  as  follows: 


In  this  way  K  must  converge  to  1.  The  normalized  vector  length  of  stress  intensity 
factors,  computed  by  the  energy  release  rate  method,  the  contour  integral  method 
and  the  cutoff  function  method,  is  shown  in  Table  5.3.  The  three  methods  are 
seen  to  converge  strongly.  Again,  the  performance  of  the  cutoff  function  method 
is  seen  to  be  the  strongest. 

Table  5.3  Normalized  vector  length  of  stress  intensity  factors  (ic) 
computed  by  the  energy  release  rate  method  (ERM),  the  contour  integral 
method  (CIM)  and  the  cutoff  function  method  (CFM) . 


p 

ERM 

CIM 

CFM 

1 

1.06576 

1.24288 

0.95868 

2 

1.02237 

0.92678 

1.03382 

3 

0.99926 

1.04100 

1.00457 

4 

0.99663 

0.98421 

1.00080 

5 

0.99658 

1.00673 

0.99967 

6 

0.99733 

0.99644 

0.99961 

7 

0.99795 

1.00104 

0.99966 

8 

0.99840 

0.99912 

0.99973 

oo 

1.00000 

1.00000 

1.00000 

5.S.  Edge  cracked  panel  problem  2. 

The  two  problems  jxist  discussed  were  constructed  so  that  only  the  first  sym¬ 
metric  and  antisymmetric  terms  of  the  asymptotic  expzmsions  were  nonzero.  This 
permitted  us  to  examine  the  performance  of  the  contour  integral  and  cutoff  func¬ 
tion  methods.  In  practical  problems  no  such  restrictions  apply  and  the  exact 
solution  is  not  known.  The  following  test  problem  is  more  nearly  representative 
of  practical  problems. 

The  problem  definition  and  mesh  desigpi  are  shown  in  Fig.  5.5.  Plane  stress 
condition  and  Poisson’s  ratio  of  0.3  are  assumed.  Solutions  were  obtained  by 


Fig.  5.5.  Edge  cracked  panel  problem  2.  Mesh  design, 
the  virtual  crack  extension  method,  the  contour  integrad  method  amd  the  cutoff 
function  method.  We  normalize  the  stress  intensity  factors  if;,  Kn  and  their  vector 
length  K  as  follows: 


&  def  Kt 

ay2irw 


V2jr  w’ 


drf  +  ■l^f 


The  results  are  shown  in  Tables  5.4  amd  5.5.  This  problem  was  solved  by  Sha  and 
Yang  using  the  virtual  crack  extension  method  [21]  and  by  Andersson  [22]  using 
a  technique  similar  to  that  proposed  by  Rybicki  and  Kanninen  [23].  The  results 
presented  herein  agree  very  closely  with  the  results  given  in  [21,22]. 

The  results  indicate  that  the  computed  values  converge  strongly.  Our  ability 
to  go  beyond  the  level  of  precision  we  actually  need  at  a  small  marginal  cost  is 
very  important  from  the  point  of  quadity  control:  There  au-e  no  other  means  for 
ensxiring  the  quadity  of  computed  data  in  practical  computations  where  the  exact 
solution  is  generally  not  known. 

Finally  we  note  that  the  contoxir  integral  and  cutoff  function  methods  require 
substamtially  fewer  CPU  cycles  thaui  the  virtual  crack  extension  method.  In  the 
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Table  5.4.  Solution  by  the  energy  release  rate  method. 
Normalized  strain  energy  and  vector  length  of  stress  intensity  factors. 


p 

N 

UIufs)E 

w 

k 

1 

43 

1.4849947 

1.439 

2 

125 

1.6162227 

1.605 

3 

221 

1.6593481 

1.587 

4 

355 

1.6940543 

1.656 

5 

527 

1.7010749 

1.666 

6 

737 

1.7030296 

1.671 

7 

985 

1.7035357 

1.673 

8 

1271 

1.7037223 

1.674 

Table  5.5.  Solutions  by  the  contour  integral  and  cutoff  fimction  methods. 
Normalized  stress  intensity  factors  aind  their  vector  length. 


P 

N 

ki 

(CIM) 

ki 

(CFM) 

kn 

(CIM) 

kii 

(CFM) 

k 

(CIM) 

k 

(CFM) 

1 

43 

0.54127 

0.42259 

-0.37480 

-0.29005 

1.650 

1.285 

2 

125 

0.49708 

0.55588 

-0.25578 

-0.28292 

1.401 

1.563 

3 

221 

0.58909 

0.56161 

-0.28951 

-0.27074 

1.645 

1.563 

4 

355 

0.57864 

0.59232 

-0.28319 

-0.29022 

1.615 

1.653 

5 

527 

0.60558 

0.59825 

-0.29398 

-0.29012 

1.687 

1.667 

6 

737 

0.59672 

0.60043 

-0.28897 

-0.29097 

1.662 

1.672 

7 

985 

0.60313 

0.60119 

-0.29196 

-0.29091 

1.680 

1.674 

8 

1271 

0.60032 

0.60132 

-0.29042 

-0.29095 

1.672 

1.674 

case  of  this  example  five  elements  have  vertices  on  the  crack  tip.  Therefore,  using 
the  central  difference  formula  for  evaluting  the  rate  of  change  of  the  potential 
energy  with  respect  to  the  crack  length,  we  needed  to  recompute  the  stiffness 
matrices  of  five  elements  twice.  At  p  =  8  the  virtual  crack  extension  method 
required  13.5  (resp.  10.5)  times  the  CPU  time  required  by  the  contour  integral 
(resp.  cutoff  function)  method.  Also,  the  contour  integral  and  cutoff  function 
methods  yield  the  Mode  1  and  Mode  2  stress  intensity  factors  separately  whereas 
the  virtual  crack  extension  method  does  not. 
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6.  SUMMARY  AND  CONCLUSIONS. 


We  have  described  two  methods  for  the  computation  of  the  amplitudes  of 
the  first  symmetric  and  antisymmetric  terms  of  the  asymptotic  expansions  of  the 
solution  at  reentrant  comers  with  stress  free  boundaries.  This  class  of  problems  is 
very  important  &om  the  practical  point  of  view  because  it  includes  the  problems 
of  linear  elastic  fracture  mechanics  and  in  many  cases  the  sites  of  failure  initiation 
2U'e  reentrant  comers. 

The  limiting  asstimptions  of  this  paper  were  adopted  in  order  to  keep  the  pre¬ 
sentation  as  simple  as  possible:  The  extraction  methods  described  herein  are  not 
intrinsically  limited  to  the  class  of  problems  considered  here,  but  implementation 
for  the  general  case  is  somewhat  more  difficult.  For  example,  we  have  considered 
real  and  simple  roots  of  (2.9a,b)  only.  If  the  solid  angle  is  360“  then  all  roots  are 
real  and  simple  and  the  procedures  presented  herein  can  be  applied  directly  to 
computing  any  number  of  coefficients  of  the  asymptotic  expansion.  If  the  solid 
angle  is  arbitrary  and  the  amplitudes  of  several  terms  of  the  asymptotic  expan¬ 
sion  are  of  interest  then  the  implementation  must  account  for  multiple  real  roots 
and  complex  roots  as  well.  Complications  arise  also  when  the  material  properties 
change,  such  as  along  the  edge  of  a  composite  panel.  In  such  cases  determination 
of  the  eigenvalues  and  eigenfunctions  is  more  difficult  than  in  the  case  discussed 
here.  (See,  for  example,  [24].)  These  difficulties  notwithstanding,  implementation 
of  procedures  for  the  computation  of  the  asymptotic  expansion  of  stress  singu¬ 
lar  terms  is  feasible.  We  have  demonstrated  that  the  amplitudes  computed  by 
the  contour  integral  and  cutoff  ftmction  methods  converge  to  their  true  values  at 
about  the  same  rate  as  the  strain  energy,  or  faster.  Thus  the  amplitudes  can  be 
computed  accurately  and  inexpensively  by  these  methods. 

The  accuracy  depends  on  the  design  of  the  finite  element  mesh  and  the  choice 
of  polynomial  degree.  In  general  the  mesh  should  be  designed  and  the  polynomial 
degree  chosen  so  that  the  desired  level  of  accuracy  is  reached  just  before  transition 
occurs  from  Phase  1  to  Phase  2  in  the  p-extension  process.  Procedures  for  correct 
mesh  design  and  proper  selection  of  p  are  discussed  in  [7,25]. 


P-extensions  can  be  implemented,  and  in  fact  have  been  implemented  in 
PROBE,  so  that  once  the  solution  for  p  =  po  is  available  the  solutions  for  p  <  po 
can  be  obtained  very  inexpensively.  This  is  because  stiffness  matrices  and  load 
vectors  have  hierarchic  structure  [17].  Consequently  the  marginal  cost  associated 
with  obtaining  stress  intensity  factors  for  p  <  ;>o  >8  very  small.  This  is  not  the  case 
when  the  h-version  of  the  finite  element  method  is  used. 

Development  of  theories  for  the  prediction  of  failure  in  metallic  and  nonmetal- 
lic  materials  is  a  problem  of  very  obvious  practical  importance.  Several  failure 
theories  have  been  proposed,  each  requiring  computation  of  some  parameters  of 
the  elastic  stress  field  (see,  for  example  [26-31]).  In  linear  elastic  fracture  mechan¬ 
ics  the  amplitude  of  the  first  symmetric  term  of  the  asymptotic  expansion  of  the 
solution  about  the  crack  tip  have  been  correlated  with  crack  extension  through 
laboratory  experiments.  It  is  conceivable  that  much  like  crack  extension,  failure 
initiation  (e.g.  the  formation  of  cracks  at  weldments  and  reentrant  comers;  the 
onset  of  delamination  in  composite  materials,  etc.)  can  be  correlated  with  the 
amplitudes  of  the  terms  of  asymptotic  expansions  also.  In  fact,  linear  elasticity 
cannot  be  useful  for  predicting  failure  initiation  events  unless  parameters  of  the 
elastic  solution  can  be  consistently  correlated  with  occurrences  of  such  events. 
Certain  parameters  can  be  computed  only  by  extraction  methods.  Others  can  be 
computed  directly  from  the  finite  element  solution  [32].  We  have  shown  that  ac¬ 
curate  and  inexpensive  determination  of  the  amplitudes  of  asymptotic  expansions 
in  the  neighborhood  of  reentrant  comers  is  possible. 
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